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ABSTRACT 

In this work we describe our implementation of a thermodynamic energy equation into the global 
corona model of the Space Weather Modeling Framework (SWMF), and its development into the new 
Lower Corona (LC) model. This work includes the integration of the additional energy transport terms 
of coronal heating, electron heat conduction, and optically thin radiative cooling into the governing 
magnetohydrodynamic (MHD) energy equation. We examine two different boundary conditions using 
this model; one set in the upper transition region (the Radiative Energy Balance model), as well as 
a uniform chromospheric condition where the transition region can be modeled in its entirety. Via 
observation synthesis from model results and the subsequent comparison to full sun extreme ultraviolet 
(EUV) and soft X-Ray observations of Carrington Rotation (CR) 1913 centered on Aug 27, 1996, we 
demonstrate the need for these additional considerations when using global MHD models to describe 
the unique conditions in the low corona. Through multiple simulations we examine ability of the 
LC model to asses and discriminate between coronal heating models, and find that a relative simple 
empirical heating model is adequate in reproducing structures observed in the low corona. We show 
that the interplay between coronal heating and electron heat conduction provides significant feedback 
onto the 3D magnetic topology in the low corona as compared to a potential field extrapolation, 
and that this feedback is largely dependent on the amount of mechanical energy introduced into the 
corona. 

Subject headings: MHD - Sun: corona 



1. INTRODUCTION AND BACKGROUND 

While the inherent complexity of the solar corona and 
its dynamic processes has engendered great scientific 
interest and study over the past half century, it is this 
complexity that is perhaps the greatest challenge when 
attempting to model and understand the underlying pro- 
cesses of observed phenomena in a quantitative fashion 
(Aschwanden 2008). For this reason, the development 
and use of fully 3D magnetohydrodynamic (MHD) 
models describing the corona and solar wind, which 
has undergone a veritable renaissance in the computing 
age, has been critical to furthering our understanding of 
the solar wind and eruptive events (Roussev & Sokolov 
2006). A particularly exciting avenue has been the 
development of 3D observation-based techniques, which 
provide the opportunity to include realistic conditions of 
the corona into studies of the global solar wind (Mikic 
et al. 1999; Roussev et al. 2003; Cohen et al. 2007) 
and to address how this governs dynamic events (Lugaz 
et al. 2007; Roussev et al. 2007). 

With its origin in the long-studied fundamental 
question of how the solar corona can be heated to and 
maintain temperatures in the million Kelvin (MK) 
regime, a primary focus of this study is to address the 
issue of energy input into the corona in the context of a 
fully 3D global model. Because non-MHD processes of 
energy transport and local conditions in the low corona 
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are responsible for coronal heating mechanisms, it is 
often quite difficult and computationally expensive to 
include the small-scale micro physics of reconnection 
and turbulence responsible for coronal heating in 3D 
models, especially since the exact processes involved 
are not universally agreed upon. As a result, heating 
models are often parameterized as a heating term that 
depends on various local magnetic and thermodynamic 
properties that is included in the energy equation (e.g. 
Aschwanden & Schrijver (2002); Schrijver et al. (2004); 
Abbett (2007); Mok et al. (2008)). However, due to 
the large number of heating models and their ad-hoc 
formulations, it is crucial to use real observations to 
compare and constrain these models and deduce the 
best combination of models, both for active regions 
and the quiet Sun, that yields the most realistic global 
simulations. When using such models to study transient 
events in the Solar Corona, particularly Coronal Mass 
Ejections (CMEs) and coronal waves seen in the extreme 
ultraviolet regime (EUV waves, Moses et al. (1997)), 
it is a motivating principle that only by constraining 
physical theories and scenarios through as many ob- 
servable manifestations available can we provide the 
best possible avenue to further our knowledge and 
precision in understanding these dynamic events (e.g. 
Lugaz et al. (2008, 2009); Cohen et al. (2009)). The 
ultimate goal being to strip away layers of empirical 
approximation with more fundamental physical terms 
and mechanisms. However, it is equally important to 
validate and asses the ability of these improvements to 
adequately represent the global conditions in the solar 
corona. 
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To address this issue through a deterministic, data- 
driven approach, we modify the global corona model of 
the Space Weather Modeling Framework (SWMF) (Toth 
et al. 2005) to include the transition region between 
the chromosphere and the corona, where non-MHD 
thermodynamic terms of energy transport, such as 
electron heat conduction, radiative losses, and coronal 
heating all become important. Using techniques similar 
to the pioneering work of Lionello et al. (2001, 2009), in 
this refined version of the model the inner boundary is 
placed at the chromosphere or upper transition region 
rather than in the low corona and the now relevant 
non-MHD terms are added to the governing MHD 
energy equation. 

In section 2 we discuss the new Lower Corona (LC) 
model and physical considerations added as part of this 
work. In section 3 we overview our runs and results and 
discuss the application of this model to the study of the 
corona (in particular the comparison to EUV imaging 
observations). In section 4 we provide a detailed analysis 
of a particular run, and benchmark it via comparison to 
a previous corona model. We conclude in section 5. 

2. THE SIMULATION TOOL 
2.1. SWMF 

The main model and starting point for this work is the 
SWMF Solar Corona (SC) model (Cohen et al. 2007). 
As far as numerical methods, the code is fully paral- 
lelized and designed to be highly customizable in terms 
the methods, solvers, and equations used. A complete 
description of the MHD equations and their implemen- 
tation in the SWMF can be found in Powell (1999) as well 
as the addition of source terms in (Groth et al. 2000). For 
the current version of the Solar Corona Model, a steady 
state solar wind solution is obtained using a variable 7 
model (Roussev et al. 2003; Cohen et al. 2007, 2008). 
The latter references use the Wang-Sheeley-Arge (WSA) 
model (Arge et al. 2004) as an initial condition to derive 
the variation in the polytropic index. This then achieves 
a steady state MHD solution after nominal integration 
(henceforth refered to as the "Standard SC model"). 

Physically, the most important advantage of this 
tool and others like it lies in its ability to simulate 
the complete 3D environment of any event and Car- 
rington Rotation (standardized solar rotation number, 
abbreviated CR) for which data for the photospheric 
magnetic field is available throughout the entire solar 
surface. The initial magnetic configuration is ex- 
trapolated using the Potential Field Source Surface 
Method (PFSSM) (Altschuler et al. 1977), which uses 
magnetic coefficients derived from observations of the 
synoptic photospheric magnetic field (typically high 
order maps using data from the MDI instrument aboard 
the SOHO observatory, or low order maps the Wilcox 
Solar Observatory). Using this method, the subse- 
quent evolution of the magnetic field towards a steady 
state during the simulation is no longer strictly potential. 

While the standard version of the SC model has been 
very successful at reproducing the bi-modal solar wind 
structure and observations at 1AU (Cohen et al. 2008), 



which was the fundamental goal of that empirical model, 
the lower boundary at the solar surface is quite smooth 
and nearly uniform. In order to reproduce and study 
the fine structures of the the low corona we significantly 
modify the SC model to address the unique physics that 
take place at this boundary (henceforth refered to as the 
the Lower Corona (LC) model). 

2.2. Including Additional Thermodynamic Terms 

This modification of the energy equation takes the form 
of: 

dE -> -> 

-QjT + V • (FmHD + F c ) = QmHD +Qr + Qh (1) 

where, the subscript MHD refers to the standard MHD 
terms in the SWMF and the additional terms F c , Q r , 
Qh are described below. Because these terms account 
for realistic energy transport, the polytropic index is no 
longer variable, and is set to a uniform value of 7 = 5/3. 

2.2.1. Heat Conduction 

The standard form of anisotropic electron heat 
conduction in the collisional Jimit is included as an 
additional enersv flux term, F c = -k T 5 / 2 B(B • VT), 
in the energy equation (equivalent physically to adding 

a source term of Q c = V • F c ). The parallel component 
of Spitzer conductivity (Spitzer 1965), this term is 
especially important in the transition region and low 
corona and plays a critical role in determining the 
equilibrium density and temperatures at the base of the 
corona. We use a value of ^0 = 1.23 x 10 -6 in cgs units. 

We also include a method for broadening the equi- 
librium scale length of the transition region, which 
is needed to resolve the entire transition region when 
using a chromospheric boundary (described in detail by 
Abbett (2007); Lionello et al. (2009)). This reduces to 

both fixing k T 5 / 2 to a constant value of fto^od anc ^ 
reducing the radiative cooling term by the inverse factor 

(T/T^ d ) for temperatures below a transition region 
value, T < T mo( i. In this work we use T mo( i = 300, 000K. 

2.2.2. Radiative Losses 

The next term, Q r = — n e n p A(T), accounts for the ra- 
diative losses of the hot coronal plasma in the optically 
thin limit. The loss function A(T) is calculated using the 
CHIANTI version 5 radiative loss routines (Landi et al. 
2006) and linearly interpolated as a tabulated function 
in the model. The assumption of a fully ionized hydro- 
gen gas is already implicit in this model, which gives 
n e n p ~ n\. It is important to note that the choice of 
abundances used in calculating the radiative losses func- 
tion can change the total losses significantly in the coro- 
nal temperature regime. For this reason we choose to use 
the coronal abundances file: sun_coronal_ext . abund 
(Landi et al. 2002) when calculating the radiative cooling 
curve with CHIANTI, and do not vary this function in 
this work. 
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2.2.3. Heating Model 1 - \B\ weighted 

The first empirical coronal heating term, Qh, that we 
examine is equivalent in form to that used by Abbett 
(2007) and was originally based on work relating the sur- 
face unsigned magnetic flux of a convective envelope to 
the x-ray luminosity of stars (Bercik et al. 2005). This 
model is expressed as: 



C(j) a 1p 

UWdV' 



(2) 



Where c = 0.8940, a = 1.1488 are fixed parameters, <p 
represents the total unsigned magnetic flux at the solar 
surface, ip is the local heating weighting function, and £ 
is a normalization constant. This formulation calculates 
the total amount of energy input into the corona and 
distributes it via -0, which is chosen to be a function of 
\B\. The primary benefit of this model is its simplicity, 
the amount of heating at any location is given by the 
value of ip, and the total amount of heat input to the 
corona (nearly constant) is normalized by the global 
integral in the denominator, J ip dV which is carried 
out over the entire domain. In this implementation, cf) is 
calculated by integrating the unsigned flux at the model 
boundary, <j> = § R=Rs \B r \dS, which is essentially the 
photospheric magnetogram used to initiate the model. 

Because both the relationship between X-Ray Lu- 
minosity and total power input into the corona is not 
well constrained, and varying magnetogram resolution 
does not constrain the true unsigned flux, ( represents a 
relatively free parameter with which to adjust the total 
power on a case by case basis. In our model runs we 
choose £ = 1/50 and linear weighting with the magnetic 
field magnitude, ip(\B\) = 1^1- This gives the heating 
term the operational equivalent of Qh = H\B\, where 
H is typically around ~ 4 x 10 -5 in cgs units, and 
total integrated power of Q tot = 3.11 x 10 28 ergs s _1 if 
applied uniformly. 

One drawback of this term (and of other \B\ weighted 
heating terms in global models) is that the heating scale 
height, defined as the point along a given loop where 
Qh(s) — Qh(s — 0)/ e ? is n °t a free parameter and varies 
with the strength and spatial distribution of B. Also, 
note that the value of H depends on the specific global 
magnetic field configuration and the radius to which the 
volume integral is calculated. In this work, we limit the 
range of influence of this function to the low corona by 
multiplying i/; by an exponential envelope function with 
a scale height of 40Mm. 

2.2.4. Heating Model 2 - exponential heating 

The second heating model is a simple exponential scale 
height model. Using a standard exponential form we 
define heating model 2 as: 



Q h = H ex V [-{r-R Q )/X}. 



(3) 



Where Hq is the local heating rate at r = R and A 
is the heating scale height. A short heating scale height 
is consistent with SOHO and TRACE EUV observations 
of coronal loops (Aschwanden et al. 2000a; Aschwanden 
& Acton 2001), and motivated by the "nanoflare" class 



of coronal heating models, in which most of the heat 
is deposited in the upper chromosphere and transition 
region (Parker 1988). It is important to note that it is 
the high thermal conductivity of the coronal plasma that 
subsequently distributes this energy along the field/flow 
lines. For heating quiet sun regions, we use the values 
Ho = 7.28 x 10" 5 ergs cm -3 s _1 , and A = 40 Mm, which 
give a total power of Q tot = 1.99 x 10 28 ergs s _1 if applied 
uniformly (surface flux = 3.27 x 10 5 ergs cm -2 s _1 ). The 
lack of explicit dependence on \B\ with an exponential 
scale height model allows us examine directly the effect 
that purely the magnetic geometry of the 3D corona has 
on the resultant thermodynamic equilibrium. 

2.2.5. Coronal hole heating 

To maintain realistic temperatures in the solar wind 
and coronal holes at large distances, we use the 
same form of eq. (3) with values of Ho = 5.0 x 
10 _T ergs cm -3 s _1 , and A = 0.7R© giving additional 
total power Q — 5.01 x 10 2T ergs s _1 . This function is 
applied uniformly in each model run. While not entirely 
realistic for the acceleration profiles observed in coronal 
hole regions and the solar wind, we find that including 
this term is adequate for reproducing coronal hole emis- 
sion in the absence of short scale height heating. 

2.2.6. Open field cutoff 

One interesting focus of this study is to isolate the 
complementary effects that non-uniform heating at 
the surface and the 3D magnetic topology have in 
determining the equilibrium structure of the corona. 
For this reason, as one of the options in the model 
we include the ability to apply the low scale height 
heating model (either 1 or 2) to only regions associated 
with closed field in the corona. To do so, we use the 
PFSSM extrapolation computed to initiate the model 
to determine at every computational element whether 
the field threading the location is open (connects 
to the source surface and becomes radial) or closed 
(turns back down to the solar surface). If the region 
is associated with open field, then only the long scale 
height "coronal hole" heating model is applied (sec- 
tion 2.2.5), while if the region is closed, both are applied. 

We would also like to note that any of the empirical 
heating terms discussed above (sections 2.2.3-2.2.6) may 
also be attributed to wave absorption and dissapation, 
probably due to the ion interaction with Alfven wave 
turbulence. In this case, the heating may be specified 
in terms of the wave (spectral) energy and their damp- 
ing coefficient, and the wave turbulence may contribute 
not only to the energy source but also to the momen- 
tum fluxes (see e.g. Sokolov et al. (2009) and papers 
cited therein). However, this particular study is focused 
primarily on the plasma conditions low corona, and we 
operate under the assumption that Alfven wave pressure 
will be much lower than the thermal pressure there. 
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2.3. Boundary Conditions 

Boundary conditions at the r = R surface are applied 
to fix constant values temperature and density at the face 
boundary according to either of the models described 
below. Each of these use a zero flow condition and fix 
the tangential component of the magnetic field (similar to 
those used in Roussev et al. (2004); Jacobs et al. (2009). 
For the supersonic flow at the outer boundary at r = 
24R we use a floating (zero-gradient) condition. 

2.3.1. Chromospheric Boundary 

The primary boundary condition used in this model 
is a simple uniform "chromospheric" condition. As de- 
scribed in Lionello et al. (2001, 2009), by setting chromo- 
spheric values of electron temperature, T e = 2 x 10 4 K, 
and density, n e = 1 x 10 12 cm -3 , and including the 
transition region broadening method described in section 
2.2.1, it becomes possible to resolve the entire transition 
region everywhere in a global model. This is a criti- 
cal advantage because the resultant topological equilib- 
rium in the coronal part of the solution (density enhance- 
ments, varying width of the transition region, etc.) will 
depend entirely on the included physics and magnetic 
geometry, and not on its proximity to the boundary of 
the model. This fact is particularly important for the 
study of any sort of dynamics in the global corona, as 
the sharp transition from fast Alfven speed in the quiet 
corona (V a ~ 2 — 3x 10 2 km/s) to extremely slow speeds in 
the chromosphere (V a ~< 50km/s), provides a realistic 
dispersion buffer for fast magnetosonic waves. 

2.3.2. REB Model: 

In order to include a realistic lower boundary condition 
for the corona, while at the same time maintaining de- 
cent computational efficiency for the global domain, we 
use the method of the Radiative Energy Balance Model 
(REB) outlined by Withbroe (1988) (further formalism, 
including coronal heating in Lionello et al. (2001)). In 
short, this model fixes the base electron temperature to 
a high transition region value, T l J = 5 x 10 5 K, and in 
turn assumes an equilibrium balance throughout the rest 
of the Transition Region below to the top of the Chro- 
mosphere, from which an electron density, nf, may be 
obtained. By rearranging the terms in the MHD energy 
equation and integrating over T at constant pressure, one 
obtains: 



\ Ko Tj{B ■ VT e ) 2 + f Q h Tj 
\\ QlThk{T)dT 



(4) 



where, the integral of the radiative loss function is car- 
ried out between Tf 1 = 10 4 and T l e r = 5 x 10 5 K. This 
allows for the base density to be a function of the 
physical conditions included in the model (coronal heat- 
ing, anisotropic heat conduction, and radiative cooling) 
and thus varies spatially along the boundary. Most im- 
portantly, the numerator of (4) shows how the interplay 

between magnetic field strength, \B\ (via Qh in the right 
term), and normal orientation at the surface, B r (via 
(B • VT e ) 2 on the left), sets the base density on the so- 
lar surface. The ability to capture these features with 



such a simple, easily computed boundary is a primary 
motivation in using this particular model. 

2.4. Geometric Considerations 

Another aspect of this work has been to address the 
unique geometric concerns that arise when including the 
transition region and below. A typical solar corona (SC) 
simulation is carried out with an adaptive Cartesian grid 
over a Sun-centered 48 x 48 x 48R cube with 4x4x4 
cell blocks. The average cell size is smallest at the 
surface (~ 10 -2 R ~ 7, 000 km) and is incrementally 
increased with distance to a largest size of 0.65 R© 
near the outer boundary, giving a few million cells. 
While this is a relatively large surface cell size, it is 
adequate considering the precision of the magnetic field 
boundary (an observational limitation). However, if one 
is required to resolve the transition region, with dynamic 
radial spatial scales on the order of a few km, it quickly 
becomes infeasible to use a Cartesian grid (with every 
1/2 refinement in length the number of surface cells 
increases by 4). This can be partially alleviated using 
methods to widen the extent of the model transition 
region to typical scales of order 300 km without affecting 
the coronal solution (section 2.2.1), but the problem is 
still quite substantial. This then represents a balance 
between height accuracy and computational efficiency 
(with a Cartesian grid, the number of surface cells 
increases by a factor of 4 with every successive 1/2 
refinement). 

To address this scale and resolution issue we adapted 
the generalized grid capability of the SWMF, which 
has ability to calculate MHD fluxes in arbitrary geome- 
tries, to construct a spherical (r, 0, (j)) grid function with 
highly non-uniform radial scales near the transition re- 
gion (dr = 230 km) that smoothly transitions to near 
equal face area at large r {dr = 30, 000 km at r = 5i?©). 
Because the grid function is continuous, the model main- 
tains the block- adaptive mesh and adaptive mesh refine- 
ment (AMR) capability for flexible local spatial reso- 
lution both initially and as the simulation progresses. 
For this study, focusing on the global structure of the 
corona, we use a uniform spacing in the angular direc- 
tions with dO = dcf) = 1.4° at the surface, which is coars- 
ened by a factor of two via AMR beyond 1.7 Rq for re- 
gions within ±65° and 1.2R Q outside (to avoid needlessly 
small cell sizes at the poles). The outer boundary is fixed 
at r = 24R . Degenerate cells touching the polar axes 
are treated with azimuthal averaging to avoid discontinu- 
ity across the poles. A grid comparison near the surface 
between this work and an AMR cartesian mesh with a 
similar number of total cells is shown in Figure 1. 

2.5. LOS image synthesis 

To study the EUV corona and associated dynamics in 
the context of real events, we include a proper treatment 
of integrated coronal emission in the EUV range, particu- 
larly for the purpose of comparing simulated observations 
to existing EUV observations by a current solar obser- 
vatory. This process, following the work of Mok et al. 
(2005, 2008); Lionello et al. (2009), involves two main 
steps: (i) characterization of detector response based on 
physical parameters that are calculated in the model and 
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(ii) line-of-sight (LOS) integration though the 3D data 
set to create synthesized image. 

The first step is to understand and calculate the instru- 
mental response for coronal material at a given temper- 
ature and density. For simplicity and ease to compare to 
studied observations we chose to model the three coronal 
band-passes at 171, 195, and 284A of the EIT instrument 
on the SOHO spacecraft (Delaboudiniere et al. 1995) and 
the AlMg filter of the SXT instrument (Tsuneta et al. 
1991). Because of the extremely high temperature and 
low densities, the emission line EUV and Soft X-Ray ra- 
diation of the solar corona above 0.5 MK can be treated 
as optically thin to a good approximation. As such, the 
intensity measured by each pixel of an imaging detector 
can be treated as a LOS integral through the coronal 
plasma, namely: 

R = J n ifi( T , n e)dl DN s- 1 (5) 

where R is the pixel response measured in Data Num- 
bers (DN) per second. Here I follows along the line-of- 
sight, the square of the electron density, n 2 e1 accounts for 
the amount of emitting material and fi(T,n e ) is the in- 
strumental response of filter i per unit emission measure 
as a function of both temperature and density. This is 
equivalent to integrating the resolved Differential Emis- 
sion Measure (DEM) distribution along I. To calculate 
the response functions, fi(T,n e ), we use the CHIANTI 
version 5 emission line analysis routines (Landi et al. 
2006) to generate synthetic spectra, and the EIT and 
SXT analysis routines for detector response (a part of 
the Solar Soft (SSW) framework written in IDL (Freeland 
& Handy 1998)). 

3. MODEL RUNS 

In this section we present the application of the LC 
model to CR 1913, initialized using a synoptic MDI 
magnetogram centered on Aug 27, 1996. Occuring 
during solar minimum, the sun at this time features a 
conspicuous equatorial coronal hole extending from the 
north to down past disk-center, as well a large active 
region. This then conveniently allows us to examine the 
ability of a given model to describe these three basic 
regimes: the average quiet-sun, coronal holes, and active 
regions. Additionally, this rotation has been studied in 
detail using similar methods (e.g. Mikic et al. (1999); 
Lionello et al. (2009)) and thus provides a means of 
comparison to existing work. 

Four model runs used to demonstrate the usefulness 
of the LC model are summarized in Table 1. As a 
convenient comparison of their global thermodynamic 
and topological properties, the LOS synthesis of EUV 
and soft X-Ray emission are shown in Figure 2. Each 
of the models presented use either \B\ weighted heating 
(model 1, section 2.2.3) or exponential heating (model 
2, section 2.2.4)). All four use coronal hole heating 
(section 2.2.5) applied uniformly to the domain, while all 
but run C use the open/closed field weighting (section 
2.2.6) to achieve a more realistic equilibrium in these 
regions. In order to include strong heating in Active 
Regions (AR) with high magnetic field strength, Run 



Table 1 

Model Runs. 



Run 


BC 


Heating 


Open Cutoff 


Power [ergs s 1 ] 


A 


REB 


oc|B| 


Yes 


2.88 x 10 28 


B 


Chromo 


<x\B\ 


Yes 


2.88 x 10 28 


C 


Chromo 


Exp 


No 


2.49 x 10 28 


D 


Chromo 


Exp + AR 


Yes 


2.88 x 10 28 



D also includes a modification to the heating function, 
and transitions smoothly from exponential heating to 
\B\ weighted heating (Q h = 4 x 10 _5 |I?| ergs cm _3 s _1 ) 
above 30 Gauss. 

All four model runs show similar global features in the 
EUV and Soft X-Ray when compared to observations 
but differ on the precise details. For example, all four 
runs resolve cool temperatures along the the polarity 
inversion lines near the surface (seen to the east and 
west of the extended equatorial coronal hole near the 
center of the disk). This is naturally due to field line 
orientation and loop length governing thermodynamic 
equilibrium. This commonality serves to emphasize the 
strong role that the 3D magnetic field topology of a 
given Carrington Rotation plays in determining dynamic 
equilibrium in the corona. 

In the comparison between Run A (REB model BC) 
and Run B (Chromospheric BC) we observe more 
uniform emission near the surface as well as slightly 
higher temperature (as seen in the longer extent of SXT 
emission) for the REB case. Because the boundary 
is near the top of the transition region in the REB 
model, the height of the transition region is less able to 
vary in response to magnitude of the downward heat 
conduction flux from above. Also, since the REB model 
does not include a region with chromospheric densities 
(where the radiative cooling is much higher due to the 
n 2 e dependence) a smaller fraction of the coronal heating 
power is lost in the transition region, leading to higher 
equilibrium temperatures. However, it is clear that the 
REB model, which is more suited for computational 
efficiency due to the longer equilibrium length scales at 
T = 500, 000K, can be used to adequately describe the 
ambient global corona. 

One of the most striking results is perhaps the 
comparison of Run C to the other model runs. In this 
run, the coronal heating function was applied uniformly 
in the computational domain (no dependence on or 
</>). While the coronal holes and significant AR emission 
are unsurprisingly not well preserved in the EUV, 
still much of the basic features of the low corona are 
reproduced. The ability of simple uniform heating model 
to reproduce obvious features on the surface, such as 
average quiet sun emission, and lessening near inversion 
lines, as well as the bi-modal structure between open 
and closed field regions at higher temperatures (284A 
line and soft X-Rays) speaks to the importance of both 
the 3D magnetic topology in determining the dynamic 
equilibrium of the corona and the redistribution of 
energy via a thermodynamic energy equation including 
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electron heat conduction. 

For a more realistic empirical model in Run D, 
we modify Run C by including AR heating and the 
open-field cutoff to better describe heating in active 
regions in coronal holes. The most interesting finding 
in this comparison is that by including significantly 
enhanced heating in the large active region, we observe 
a pronounced effect on the equilibrium structure of 
the AR associated closed field/streamer region (East- 
ern side of the corona in Figure 2). This produces 
noticeable feedback on the global structure as seen 
via a larger northward tilt of the heliospheric current 
sheet near the AR longitudes (shown in detail Figure 
3), and suggests a complex relationship both between 
the thermodynamics of the low corona and the global 
structure of the solar wind. Because this effect cannot 
be fundamentally extracted from extrapolating the 
magnetic field alone, this provides not only a stong 
motivation for using thermodynamic models, but also 
another element with which to constrain and refine 
AR heating models in future studies (the one used in 
in this work being a necessary but crude approximation) . 

Another important inference about coronal heating can 
be gleaned from the comparison of the two heating mod- 
els directly (Run B and Run D). Although the total heat- 
ing power is the same in each model, the \B\ weighted 
heating model naturally has a shorter scale height (due 
to the falloff of \B\ with radius) and we subsequently ob- 
serve slightly enhanced pressure in the low corona and 
lower temperatures (hence density scale height) in the 
closed field regions. However, the fact that the simple ex- 
ponential heating model applied to the quiet sun in Run 
D can reproduce the mean structure of the low corona 
with the same fidelity as a more complicated empiri- 
cal model, demonstrates that perhaps complex empirical 
models are not as entirely necessary as some might sug- 
gest and motivates a need for more physics based models 
to advance further. 

4. DETAILED ANALYSIS 

In this section we provide a more detailed analysis of 
Run D in particular. We examine the magnetic and 
thermodynamic structures realized by the LC model and 
provide a comparison to the standard SC model (Cohen 
et al. 2007). 

4.1. Temperature structure of low corona 

An important advantage of including heat conduction 
in a global 3D enrionment is the ability to the study 
complex open/closed field topologies in a self consistent 
manner, which allows one to study the feedback effect 
that the magnetic topology has on heat distribution and 
vice versa. The most striking effect can be seen by ex- 
amining the differences between heat conduction along 
open or closed field geometries. Because closed field re- 
gions represent closed systems for the flow of thermal en- 
ergy, they acheive higher temperatures than their open- 
field counterparts and shift into the soft-Xray regime 
(T ~> 1.5MK). This natural correspondence can be 
clearly seen in Figure 4 where we show a 3D surface at 
fixed temperature (T e = 1.6 MK) and an LOS image of 
soft-X-Ray emission. 



4.2. Magnetic structure of low corona 

It is also important to examine in detail the effect that 
the thermodynamic model has on the magnetic structure 
as a whole. In Figure 5 we show a comparison between 
closed field lines at r = 2. OR© for the PFSSM initial con- 
dition (Altschuler et al. 1977), the standard SC model 
(Cohen et al. 2007) and Run D. Immediately obvious 
are the drastic changes between the closed field/streamer 
structure of the polytropic SC model and the thermo- 
dynamic LC model. Unsurprisingly, the high tempera- 
tures and thermal energy introduced by coronal heating 
and conduction along the closed field produces a large 
thermal stress in the relatively high f3 regions near the 
streamer cusps along the current sheet. This in turn 
leads to significant reorganization of the magnetic field 
as compared to a PFSSM or polytropic MHD model. The 
fact that this sort of stress is highly sensitive to both the 
local nature of a given coronal heating model (Figure 3) 
and the global magnetic topology (Figure 5) is a leading 
motivation for using a fully 3D, thermodynamic model 
when studying the structure of the global corona. 

4.3. EUV and Soft X-Ray comparison 

In Figure 6 we show LOS EUV and Soft-X-Ray image 
synthesis comparisons between run D and the standard 
SC model. In Figure 7 we show two quantitative slice 
comparisons cutting an average quiet sun region and the 
large active region. The improved agreement of Run D 
in multiple filter bands is indicative of the improved tem- 
perature resolution possible with the LC model. While 
not an entirely fair comparison because the standard 
SC model does not include a full energy equation, we 
demonstrate the ability to quantitatively asses the rela- 
tive agreement between a given model and observations 
in the EUV regime. This sort of comparison also allows 
one to study in detail various choices of heating models 
and, though multi-filter comparison, provide a basis to 
resolve the inherent degeneracies between modifying the 
total coronal heating power and the heating scale height 
within the model. A further examination of some of the 
limitations imposed by EUV synthesis from MHD mod- 
els with finite resolution as well as future considerations 
can be found in Appendix A. 

4.4. Thermodynamic Comparison 

It is also instructive to examine the global temperature 
and density structure as a function of height for both 
Run D and the standard SC model. In Figures 8 and 
9 we show electron temperature and number density 
at three separate heights in the low corona. Important 
to emphasize is the fact that with a thermodynamic 
model, it is possible to resolve both significant changes 
in temperature with height and with respect to the 
magnetic field topology in the model. As expected, near 
the top of the transition region Run D resolves cooler 
temperatures near the inversion lines due the short 
loop length, while simultaneously higher in the corona, 
hotter temperatures (T > 1.5MK) are achieved in close 
field regions, a natural result of energy flowing on a 
closed path via heat conduction. This greater sensitivity 
to field geometry (something difficult for a polytropic 
model to do without extreme fine tuning of boundary 
conditions) is critical when trying to correctly describe 
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the conditions in the low corona. 



5. CONCLUSION 

In this paper we demonstrate the clear need for a 
self-consistent treatment of the thermodynamic energy 
equation and boundary conditions when studying the 
global properties of the low corona with an MHD model. 
We observe that the interplay between coronal heating 
and electron heat conduction strongly governs the 
details of important strucures in the low corona, which 
is something that connot be described by extrapolation 
of the magnetic field on its own (e.g. the significant 
stress of the streamer regions and current sheet in 
response to active region heating, Figures 3 and 5). We 
also showcase the ability of the LC model to examine 
the effectiveness and applicability of various heating 
functions, and find that a simple exponential heating 
profile does quite well in describing the average structure 
of the quiet sun. 

With this foundation, we believe that the LC model 
can provide a more realistic vehicle with which to study 
a wide range of aspects of the dynamic corona and so- 
lar wind. For example, the time-dependent dynamics 
of any sort of transient event in the corona, such as a 
CME or EUV wave, depend critically on the magnetic 
and density structure of the corona (density scale height 
being highly sensitive to temperature) and therefore it 
is important to use improved models when determining 
the global ambient medium in which they evolve. Addi- 



7 

tionally because of the general flexibility of the SWMF, 
future coupling of the LC model to complementary so- 
lar wind models, such as the SC model based on a full 
description of Alfven wave turbulence currently in de- 
velopment (Oran et al. 2009, in press), or one including 
a multi- species treatement to resolve electron/ion tem- 
perature decoupling in the high corona, has become a 
feasible scenario. Ultimately, one thing is clear: with 
each passing month, year, and decade, our understand- 
ing of these fundamental process is constantly refined. 
In furthering this process it is critical to continue to im- 
prove the existing models and validate new efforts as they 
emerge. 
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APPENDIX 

INTERPRETING SYNTHESIZED EUV IMAGES FROM GLOBAL MODELS 

One critical caveat when imaging the solar corona is that remote sensing images are fundamentally ambiguous in 
the sense that they are 2D realizations of a fully 3D, optically thin data set. For our purposes, it is best to frame 
the problem with the concepts of modeling an unresolved Differential Emission Measure (DEM) distribution and 
Hydrostatic Weighting Bias, well described in Aschwanden & Nitta (2000b); Aschwanden & Acton (2001). These 
issues arise due to two main facts, (i) The spatial resolution of current imaging instruments (e.g. TRACE, SOHO, 
STEREO) all lie below the fundamental width of an isolated flow/field line. And (ii) any image is a line of sight 
integral over a range of heights, which are preferentially biased towards hotter flow lines with longer density scale 
heights. Thus at further distance off of the limb, the observation may be biased towards flow lines that may not be 
energetically important at the surface. Local DEM distributions can also be directly observed through coupled filter 
and tomographic inversions (Frazin et al. 2009). 

Though the LOS integral for image synthesis (5) includes the range of temperatures and densities in cells along the 
line of sight, it cannot fundamentally account for the DEM distribution that is below the physical resolution of the 
model. In MHD, a given cell has a unique value for its thermodynamic properties, particularly n e , and T, while a 
physical picture of the low (3 multi-hydrostatic corona below the resolution limit can be treated as a superposition of 
densities and temperatures, each confined to an isothermal flow/field line. In this sense, the true DEM distribution of 
a parcel of gas that contributes to an EUV image (e.g. the EIT instrument) is instead represented as a delta function 
in the simulation. One aspect of this limitation can be seen in difference between the data and synthesized limb 
profiles (Figure 7). Because the observational weighting bias shifts to higher temperatures further with distance from 
the surface, the observed limb profile cannot be represented in all four filters simultaneously with only the larger scale 
temperature equilibrium produced by the MHD model. 

To further examine this issue in the context of the LC model we use filter ratios between the EUV images, which are 
nominally sensitive to temperature differences. In Figure 10 we show filter ratio values for SOHO EIT observations on 
Aug 27, 1996 and the path in the filter ratio plane spanned by all possible filter ratios of the optically thin response 
functions fi$§(T,n e )/ fm(T,n e ) and f284(T,n e )/ fi$5(T,n e ) (section 2.5). Unsurprisingly, the EUV data spans a 
much larger range than that allowed by a single temperature. To simply examine the effect that an unresolved DEM 
distribution may have on the synthesis and subsequent comparison to observational data we construct an ad hoc DEM 
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by Gaussian convolution of the filter response functions with T: 

i r 10? (t-t ) 2 v 

fr d (T ,n e ) = ^— fi(T,n e )e * 2 t } dT (Al) 

and use gt = 0.5To for this discussion. While we by no means posit that this DEM distribution should resemble the 
corona everywhere, it includes the basic temperature splitting features of those gaussian DEM distributions calculated 
from observations (e.g. Aschwanden & Acton (2001); Frazin et al. (2009). We also show the effect of this modification 
in the right panel of Figure 10. 

The effect that this simple induced DEM distribution has on the synthesized ratios is shown in Figure 11. It is clear 
comparing the limb tracks of the synthesized images to the temperature tracks in Figure 10 that the LOS synthesis 
at the current model resolution does not produce a significant DEM distribution on its own (particularly at the limb). 
However, the improvement gained in terms of the relative location with respect to observations for the modified DEM 
synthesis suggests that this effect should not be overlooked when constructing more detailed models in the future. 
High resolution studies designed to approach the fundamental scale of isothermal separation between coronal loops, 
particularly those involving local, time-dependent heating for example, should be able to address this fundamental 
issue. 
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Figure 1. Comparison of grid geometry and cell sizes near the solar surface for the SWMF SC model (coloring indicates radial 
extent). Left: a Cartesian grid with minimum ds — 6, 800 km. Right: The customized spherical grid with minimum dr — 230 km 
(part of this work). All three directions do not have to be further refined at the surface to enhance the radial resolution when 
using spherical geometry, leading to significant savings in computational overhead. In this example, both simulation domains 
have comparable numbers of total cells. 
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Figure 2. Comparison of EIT 17lA, 195A, 284A, and SXT AlMg image synthesis to observations for the model runs of CR 
1913 centered on Aug 27, 1996. Runs A, B, C, and D are shown top to bottom. Bottom row: SOHO EIT and SXT Almg 
observations of the same date near 01:00 UT (00:10:13, 00:24:05, 01:05:19, and 01:07:28 respectively). Run A: REB BC + B 
weighted heating + Open Field Modification. Run B same as A with chromospheric BC. Run C: Chromospheric BC + uniform 
exponential heating. Run D: chromospheric BC + Exponential heating + B weighted AR component + open field modification. 
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Figure 3. Comparison of the significant effect that AR heating has on the tilt of the heliospheric current sheet between Run C 
(uniform heating, left) and Run D (with AR heating and open field modification, right). The thermodynamic stress induced by 
significant local heating of the active region (~ 7x 10 27 ergs s -1 ) causes an extreme shift in the northward tilt of the heliospheric 
current sheet ^35°. Simultaneously, the region not thermally connected to the active region on the opposite side of the sun 
remains unchanged. The color contours display the radial component of the magnetic field, B r . The chosen slice is parallel 
to the polar axis and intersects the inversion line between the large active region seen on Aug 27, 1996. The 3D magnetic 
streamlines intersect the slice at the same location in each figure. 
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Figure 4. Visual comparison of 3D topological structure (left) to LOS synthesis of SXT AlMg response (right) for Run D 
(Exponential heating + AR heating). The surfaces shown are of constant r = 1.03 R© and T = 1.6 MK with color contours of 

T. The black contour lines of normal radial field direction \B r \ — 0.2 are shown to illustrate how the hot, closed field streamer 
regions overlay the inversion lines on the surface. 
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Figure 5. Magnetic field comparison of the streamer structure of CR 1913 between 3 different models, from a top down (top) 
and Earth-centered viewpoint (bottom) at 01:00 UT Aug 27 1996. Left: The PFSSM initial condition (Altschuler et al. 1977). 
Middle: Steady state equilibrium of the default SC model (Cohen et al. 2007). Right: Steady state of Run D (Exponential 
heating + AR heating). It is immediately clear that the closed field structures become highly stressed by heating near the 
surface and the subsequent redistribution via heat conduction. Selected field lines are chosen at near equal intervals along the 
intersection of r = 2. OR© and \B r \ = 0.5 in each model. 
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Figure 6. Comparison of EIT 17lA, 195A, 284A, and SXT AlMg image synthesis to observations for Aug 27, 1996. Top 
row: steady state of the standard SC model (Cohen et al. 2007). Middle row: Run D (Exponential heating + AR heating), 
demonstrating significantly improved agreement with surface conditions and streamer topology. Bottom row: SOHO EIT and 
SXT Almg observations of the same date near 01:00 UT (00:10:13, 00:24:05, 01:05:19, and 01:07:28 respectively). 
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Figure 7. Quantitative slice comparison of EIT 171 A, 195A, 284A, and SXT AlMg synthesis from Run D (blue), the standard 
SC model (green), and observations (orange). The slices are chosen to include the average quiet Sun (left) and a large active 
region (right) observed on the disk. 
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Figure 8. Temperature contours at r = 1.01 (top), 1.03 (middle), and 1.10 (bottom) R© for run D (left) and the default SC 

model (right) for CR 1913 (Black lines at \B r \ = 0.2). The full thermodynamic energy equation is able to capture large changes 
in temperature over short changes in radius. As expected, near the top of the transition region we resolve cooler temperatures 
near the inversion lines due the short loop length. Simultaneously higher in the corona, hotter temperatures (T > 1.5MK) are 
achieved in close field regions, a natural result of energy flowing on a closed path via heat conduction. 
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Figure 9. Electron density contours at r = 1.01 (top), 1.03 (middle), and 1.10 (bottom) R0 for run D (left) and the standard 

SC model (right) for CR 1913 (Black lines at \B r \ — 0.2). This demonstrates the high dynamic range in density over small 
changes in radius resolved with the LC model. 
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Figure 10. EIT filter ratio path with no DEM modification (left) and ratio path with Gaussian convolution of response with 
(tt = 0.5Tb (right). The blue and green points represent EUV data on an off the limb (binned to 1000 points each). The 
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curve path follows the ratio values [/i9s(T, n e )//i7i (T, n e ) and f284(T, n e ) / f 195 (T,n e )] from T = 2.5 x 10 
of the curves represents the variation of the ratio values from n e = lx 10 6-12 cm -3 at a given temperature. This is intended to 
illustrate the large variation of filter ratio values with respect to an unresolved DEM distribution. 
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Figure 11. Effect of simple DEM modification on the synthesized EIT filter ratios for Run D. LOS synthesis ratios with no 
DEM modification (left) and those with Gaussian convolution of response with ar — 0.5To (right). While the location of the 
synthesized filter ratios improves with this DEM modification, the fine details of the local DEM distribution, realized by the 
variation of observed EIT filter ratio values, are difficult to reproduce with an MHD model involving large, global scales. 



